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A GAS-KINETIC BGK SCHEME FOR THE COMPRESSIBLE NAVIER-STOKES 

EQUATIONS* 

KUN XUt 


Abstract. This paper presents an improved gas-kinetic scheme based on the Bhatnagar-Gross-Krook 
(BGK) model for the compressible Navier-Stokes equations. The current method extends the previous gas- 
kinetic Navier-Stokes solver developed by Xu and Prendergast by implementing a general nonequilibrium 
state to represent the gas distribution function at the beginning of each time step. As a result, the requirement 
in the previous scheme, such as the particle collision time being less than the time step for the validity of 
the BGK Navier-Stokes solution, is removed. Therefore, the applicable regime of the current method is 
much enlarged and the Navier-Stokes solution can be obtained accurately regardless of the ratio between 
the collision time and the time step. The gas-kinetic Navier-Stokes solver developed by Chou and Baganoff 
is the limiting case of the current method, and it is valid only under such a limiting condition. Also, in 
this paper, the appropriate implementation of boundary condition for the kinetic scheme, different kinetic 
limiting cases, and the Prandtl number fix are presented. The connection among artificial dissipative central 
schemes, Godunov-type schemes, and the gas-kinetic BGK method is discussed. Many numerical tests are 
included to validate the current method. 

Key words, gas-kinetic method, Navier-Stokes equations, Chapman-Enskog expansion, kinetic bound- 
ary condition, artificial dissipation, Godunov method 

Subject classification. Applied Numerical Mathematics 

1. Introduction. There are many approaches for the numerical solution of the compressible Navier- 
Stokes equations. Godunov-type schemes solve the Navier-Stokes equations in two steps, i.e., the inviscid 
Euler step and the viscous step. The Euler solution is based on an exact or approximate Riemann solvers. 
For the viscous part, a central difference method is generally adapted [13, 31]. 

Based on the gas-kinetic theory, the Navier-Stokes equations can be derived from the Boltzmann equation 
using the Chapman-Enskog expansion. Therefore, a Navier-Stokes solver can be equally obtained by solving 
the Boltzmann equation, especially the simplified collision models [3, 2]. In the gas-kinetic representation, 
all flow variables are moments of a single particle distribution function. Since a gas distribution function 
is used to describe both equilibrium and nonequilibrium states, the inviscid and viscous fluxes are obtained 
simultaneously. Furthermore, due to the Boltzmann equation, the kinetic method and the Direct Simulation 
Monte Carlo (DSMC) method could possibly be matched in the near continuum regime [18]. But, this 
does not mean that the gas-kinetic schemes are always superior in comparison with Godunov-type schemes 
for the Navier-Stokes solutions. There is also an operator splitting procedure in solving the Boltzmann 
equation. In many kinetic schemes, the free transport equation or the collisionless Boltzmann equation, i.e., 
ft+uf x = 0, is used for the flux evaluation across a cell interface. Then, the collision part, i.e., ft = Q(f, /), 
is implemented inside each cell. Even though a nonequilibrium gas distribution function / 0 can be used 
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flCASE, Mail Stop 132C, NASA Langley Research Center, Hampton, VA 23681-2199 (email: kxu@icase.edu), and Mathe- 
matics Department, the Hong Kong University of Science and Technology (email:makxu@uxmail.ust.hk). 
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as the initial condition for the free transport equation [6], as discussed in the current paper the validity of 
this kind of kinetic methods for the Navier-Stokes equations requires p/p At, where p is the dynamical 
viscosity coefficient, p is the pressure, and At is the numerical time step. With the introduction of particle 
collision time r, the above requirement is equivalent tor> At. Note that for any viscosity coefficient p, the 
corresponding particle collision time r can be obtained from a simple kinetic model [22, 17]. The underlying 
reason for the above requirement is that the free transport mechanism introduces a numerical collision time 
At, which generates a numerical viscosity being proportional to it, i.e., p n ~ At [19]. Therefore, the condition 
t At means that the physical viscosity coefficient (~ r) should be much larger than the numerical one 
(~ At) in order to have an accurate Navier-Stokes solution. As analyzed in this paper, in some situations 
the above requirement cannot be satisfied, and the above kinetic scheme can perform poorly. 

The BGK scheme differs from the above kinetic method is mainly on the inclusion of particle collision 
time t in the gas evolution stage. Instead of solving the collisionless Boltzmann equation, a collisional BGK 
model is solved for the flux evaluation, i.e., ft+uf x = (g — f)/r [1]. As a consequence, the dissipation in the 
transport process is controlled by the collision time r instead of the time step At. As analyzed previously 
[42, 37], the BGK scheme does give the Navier-Stokes solution in the region where r < At. Under such a 
situation, the distribution function used for the flux evaluation in the BGK method automatically goes to 
a Chapman-Enskog expansion of the BGK equation. The current paper is about the extension of the BGK 
scheme by including a non-equilibrium state as the initial condition of the gas distribution function. As 
a result, a nearly consistent kinetic method for the Navier-Stokes equations is developed, which is valid in 
both t < At and t > At cases. In other words, the influence of time step At on the accuracy of the viscous 
solution is reduced to a minimal level. Many test cases are included to support the arguments. Also, in the 
current paper, the appropriate implementation of boundary condition, Prandtl number fix, and the relation 
among the schemes with arificial dissipation, upwinding, and kinetic approximation, will be discussed. 

2. A BGK Scheme. The fundamental task in the construction of a finite-volume gas-kinetic scheme 
for the compressible flow simulation is to evaluate a. time-dependent gas distribution function / at a cell 
interface, from which the numerical flux can be obtained. 

2.1. Reconstruction. Following van Leer’s MUSCL idea [32], a numerical scheme is composed of 
an initial reconstruction stage followed by a dynamical evolution stage. At the beginning of each time step 
t = 0, cell averaged mass, momentum and energy densities are given. For a higher order scheme, interpolation 
techniques must be used to construct the subcell structure. Simple polynomials usually generate spurious 
oscillations if large gradients exist in the data. The most successful interpolation techniques known so far 
are based either on the TVD, ENO or LED principles[10, 11, 15]. These interpolation techniques can be 
applied to the conservative, characteristic or primitive flow variables. In this paper, the reconstruction is 
solely applied to the conservative variables. The limiter used is the van Leer limiter. With the cell averaged 
conservative variables Wj, and their differences s + = {wj + 1 — wj)/Ax and s_ = (wj — Wj- i)/ Ax, the slope 
of w in cell j is 

L(s + ,s.) = S(s + ,s _) |S +^- 1 
|s+| + I s-\ 

where 5(s + ,s_) = sign(s + ) + sign(s_). After reconstruction, the conservative variable w inside cell j is 
distributed linearly, 

wj(x) = wj + L(s + ,s-)(x — Xj), 
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and the interpolated flow distribution around a cell interface is shown in Fig. (4.1). The BGK scheme is 
basically to present a numerical Navier-Stokes solution from the above macroscopic initial condition, where 
the inviscid and viscous fluxes are obtained simultaneously in the evolution of the gas distribution function 
/• 

2.2. BGK Model. Since we are going to use a directional splitting method to solve the 2D BGK 
equation. The BGK model in the x-direction can be written as [17] 

( 2 . 1 ) ft + ufx = ~ — — > 

T 

where / is the gas distribution function and g is the equilibrium state approached by / . Both / and g are 
functions of space x, time t , particle velocities («,«), and internal variable £. The particle collision time r is 
related to the viscosity and heat conduction coefficients. The equilibrium state is a Maxwellian distribution, 

g = p^-)^ e -K{y-uf+(v-vr+e) ^ 

v 7 r 

where p is the density, U and V are the macroscopic velocities in the x and y directions, and A is related 
to the gas temperature m/2kT. For a 2D flow, the particle motion in the 2 direction is included into the 
internal variable £, and the total number of degrees of freedom K in £ is equal to (5 — 3q)/(7 — 1) + 1. 
In the equilibrium state, £ 2 is equal to e = $+& + -+&■ The relation between mass p , momentum 
(n = pU, m = pV), and energy E densities with the distribution function / is 


( 2 . 2 ) 



a = 1,2, 3, 4, 


where ip a is the component of the vector of moments 

= W’i,V’ 2,V’3,V’4) T = (l,u,v, ^(u 2 +v 2 + £ 2 )) t , 

and dE = dttdvdf is the volume element in the phase space with df = dfi df- 2 ---dfK- Since mass, momentum 
and energy are conserved during particle collisions, / and g satisfy the conservation constraint 

(2-3) J{g- f)ip a dS = 0, a = 1, 2, 3, 4, 

at any point in space and time. For an easy reference, the formula of the moment of a Maxwellian in 2D are 
presented in Appendix A. 

For a local equilibrium state with / = g, the Euler equations can be obtained by taking the moments of 
ip a t° Eq.(2.1). This yields 


/ 


ipaigt + ug x )d E = o, 


a = 1 , 2 , 3 , 4 . 


and the corresponding Euler equations in ,r-direction are 


/ p \ 


( PU \ 

pU 

+ 

pU 2 + p 

pV 

pUV 

V E ) 

t 

\(E + p)u) 


3 



where E = \p{U 2 + V 2 + &&) and p = p/2X. 

On the other hand, to the first order of r, the Chapman-Enskog expansion gives / = g — r(gt + ug x ). 
Taking moments of ip again to the BGK equation with the new /, we get 

J ip{gt+ ug x )dE = t J ip(g tt + 2 ug xt + u 2 g xx )dZ , 

from which the Navier-Stokes equations with a dynamic viscous coefficient p = rp can be obtained, 


(2.4) 


where 


/ p \ 


( pU \ 


( 0 \ 

pU 

+ 

pU 2 +p 



pV 

pUV 


$2x 

V E ) 

t 

\ (E + p)U ) 

X 

\ $3x J 


^la: 


r dU 
tp[2 & 


2 dU <9F 
K + 2^dx + dy )J ’ 


,dV dU s 

S2x=Tpi d^ + ^ 


[nTT du jr ,dv du , 2 

® 3 * - Tp[ U q x + V ( Qx + Qy ) R + 2 


TT ,dU dv , K + 4 3,1, 
U{ dx + dy )+ 4 dx^A } 


For the ID flow, where only {/-velocity exists, the above viscous governing equations become 


( P\ 

pU 

\Ej 


( pU \ 

pU 2 +p 

\(E + P)UJ X 


( 


0 

2K rrpU x 


\ 


K + 1 

^ T P a) x + ^ 


K + 1 T PUUx ) x 


where K = (5 - 37)/ (7 - 1) and E = \p{U 2 + From the above Navier-Stokes equations, an exact 

solution of a shock structure can be obtained. 


2.3. BGK flow solver. The general solution of / of the BGK model at a cell interface Xj+ 1/2 and 
time t is 

1 ft 1 

(2.5) f(x j+1 / 2 ,t,u,v,0 = - / g{x' ,t' ,u,v,i)e~ (t ~ tl),T dt' +e~ t,T f 0 {x j+ i/ 2 -ut), 

T Jo 

where x' = Xj+ 1/2 — u(t — t ') is the trajectory of a particle motion and fo is the initial gas distribution 
function / at the beginning of each time step (t = 0). Two unknowns g and fo must be specified in Eq.(2.5) 
in order to obtain the solution /. In order to simplify the notation, aq+1/2 = 0 will be used in the following 
text. 

In all previous BGK schemes [37], based on the initial macroscopic variables, see Fig. (4.1), the initial 
gas distribution function / 0 is assumed to be 


(2.6) 


( g l [ 1 + a l x] , x < 0 
\ g r [1 + a r x\, x > 0 


where g l and g r are the Maxwellian distributions at the left and right of a cell interface. The slopes a 1 and a r 
are coming from the spatial derivative of a Maxwellian and have a unique correspondence with the slopes of 
the conservative variables. Note that the formulation of a 1 and a r will be given later. The basic assumption 
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in the above formula is that, even with a discontinuity at the cell interface, the gas is assumed to stay in an 
equilibrium state on both sides of the discontinuity. This assumption is valid for any flow simulation, where 
the cell size Ax cannot properly resolve the viscous flow structure, such as in the shock capturing case of 
the Euler equations. When the cell size is much larger than the shock thickness, the shock does appear as 
a discontinuity and the flows in the upstream and downstream stay in equilibrium states. However, if the 
mesh size is fine enough to well resolve the physical shock structure, the initial gas distribution function / 0 
should give an accurate description of the real physical situation inside a shock wave, which deviates from an 
equilibrium Maxwellian. Therefore, a non-equilibrium state must be used to represent the physical reality 
in this case! So, in order to represent a general situation, in the current paper the initial gas distribution 
function / 0 will be assumed to have the form, 

( g l [1 + a l x - r(a l u + A 1 )] , x < 0 
1 g r [1 + a r x — r(a r u + A r )], x > 0 

where additional terms represent the nonequilibrium states from the Chapman-Enskog expansion of the BGK 
model. Again, the detail formulation of (a l ,A l ,a r ,A r ) will be given at a later time. Basically, the additional 
terms of —r(a l u + A l )g l and —r(a r u + A r )g r account for the deviation of a distribution function away from 
a Maxwellian. Since the nonequilibrium parts have no direct contribution to the conservative variables, i.e., 


f (a l u + A l )i/>g l dZ = 0, 
I ( a r u + A r )*l>g l dE = 0, 


both distributions Eq.(2.6) and Eq.(2.7) represent basically the same macroscopic distributions shown in 
Fig. (4.1) ! Hence, we can clearly observe that a gas-kinetic approach does have more freedom to describe a 
flow. To keep an initial non-equilibrium state in the gas distribution is physically necessary and numerically 
possible. It gives a more realistic description of the flow motion in the dissipative region. Many kinetic 
schemes have used the Chapman-Enskog distribution function as the initial condition [6, 16]. 

After having / 0 , the equilibrium state g around (x = 0,f = 0) is assumed to have the same form as that 
proposed in the previous BGK schemes [37], 

(2.9) g = go [l + (1 — Yi.[x])d l x -I- H[a;]d r x + At] , 


where Hfarl is the Heaviside function defined as 


m = { 


x < 0, 
x > 0. 


Here go is a local Maxwellian distribution function located at x = 0. Even though, g is continuous at x = 0, 
but it has different slopes at x < 0 and x > 0, see Fig. (4.2). In both f 0 and g, a 1 , A 1 ,a T , A r ,a l ,a T , and A 
are related to the derivatives of a Maxwellian in space and time. 

The dependence of a l ,a r ,...,A on the particle velocities can be obtained from a Taylor expansion of a 
Maxwellian and have the following form, 

a 1 = a[ + a l 2 u + a l 3 v + <4^(w 2 + w 2 + £ 2 ) = o' a ^ a , 


A 1 = A[ + A , 2 u + A\v + A l 4 -(u 2 +v 2 +£ 2 ) = a l a ip a , 
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A — Ai + A 2 'U + A 3 V + A 4 ~(u 2 + £ 2 ) — A a 1 p a , 

where a = 1,2, 3, 4 and all coefficients a\ , a 2 , ..., A 4 are local constants. 

In the reconstruction stage described earlier, we have obtained the distributions pj{x), rhj(x), fij{x), 
and Ej(x) inside each cell Xj_ i/ 2 < x < Xj +1 / 2 . At the cell interface Xj +i / 2 , the left and right macroscopic 
states are 

/ Pj(x j+1 / 2 ) \ / Pj+l{Xj + i/ 2 ) \ 


Wj{x j+1/2 ) = 


m j (x j+ 1 /2 ) 

flj (Xj + i/ 2 ) 

\ Ej(xj +1 / 2 ) ) 


, Wj+i(Xj+ 1/2) — 


m j+ 1 ( x j+ 1/2) 
n J+l (x 1+ i/ 2 ) 

V A'j+l ( x j+l/2) ) 


By using the relation between the gas distribution function / and the macroscopic variables (Eq.(2.2)), 
around Xj +i / 2 we get 


(2.10) 


( 2 . 11 ) 


f i IJ ~ _ , f 1 1 (J _ Wj(*i+i/ 2 )-t«i(«i) 

/ 5 i/xfc, = Wj(*j+i/2) ; / fir a 


J g r ipdE = w j+1 {x j+1/2 ) ; J g r a' 


_ Wj+i(Xj+i) - U>J+1 (x j+ i /2 ) 


V>eG = 


Ax+ 


where Ax = Xj +1 / 2 - Xj and Ax + = Xj + 1 - Xj +1 / 2 . With the definition of the Maxwellian distributions, 

-A'((u-£p) 2 +(u-V') 2 +C 2 ) 


K+ 2 

\l 2 

9 l = P l (~) e 
7 r 


g r =p r ( — ) 2 e -A’’((w-t/’') 2 + (»'- y ’') 2 +4 2 ) ) 

7 T 


and from Eq.(2.10) and (2.11), all the parameters in g l and g r can be uniquely determined, 


and 


where 


and 


/ A 


^ Pi i x j+i/2) ^ 

u l 


ffij (.X i+ i/ 2 )/ pj (Xj_)_i/ 2 ) 

v l 


™j{ x j+ 1/2)/ Pi (^+1/2 ) 

V A ‘ J 


V A ' / 



Pi+i (^+1/2) 


U r 

v r 

\ Ar ! 


A' = 


m j + 1 ( x j+ 1/2)/ Pj + 1 ( x j+ 1/2) 
%+i (*3+1/2)/ Pi+i (^i+1/2) 

y A r 

(/f + 2)pj(x J+ i/ 2 ) 


A r = 


4 (-®j(*j+i/2) - |(m|(x j+ i /2 ) + n‘j(xj + i/ 2 ))/ pj(xj + i/ 2 )) 

(AT + 2 )pj + 1 (xj +1 / 2 ) 

4 {Ej+l{ x j + l/ 2 ) ~ \{fh 2 j +l { Xj+ 1 /2 ) + fi 2 +1 (Xj+1/2))/ Pj+l i x j + i/ 2 )) ‘ 
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Once g r is obtained from the above equations, the slope a r in Eq.(2.11) can be computed from, 


(2.12) 


(al\ 


Wj+l(£j+l) Wj+l {Xj+1/2 ) 
p r Ax+ 


= m: 


a0 


l 3 

W 


= M r a0 a 


0 ' 


where M* 0 = f g r ip a ip 0 dE/ p r . The matrix and the direct evaluation of the solution ( al,a 2 ,al,a\) T from 
the above equation are presented in Appendix B. For g l , the matrix M l a0 — J g l 4’ a 4’0dE/ p l has the same 
structure as M* 0 , (a[,a l 2 ,a l 3 ,a l 4 ) T in Eq.(2.10) can be obtained similarly using Appendix B. After having 
the terms a 1 and a r , A 1 and A r in / 0 can be found from Eq.(2.8), which are 


M a0 A 0 — J ci uip a d~, 

(2.13) M r a0 A r 0 = ^ / a r ui’ a dE. 

P J 

Since M l a0 , M* 0 , and the right hand sides of the above equations are known, all parameters in A 1 and A r 
can be obtained subsequently using the method in Appendix B again. 

After determining / 0 , the corresponding values of po,Uo,V 0 and A 0 in g 0 Eq.(2.9), 


K+ 2 

g 0 = p 0 ( -) 2 e- x °« u - u ^ + ( v - v rf + ^ 

7 r 

can be determined as follows. Taking the limit t ->■ 0 in Eq.(2.5) and substituting its solution into Eq.(2.3), 
the conservation constraint at (a: = Xj + i/ 2 ,t = 0) gives 

(2.14) [g 0 ipdE = w 0 = [ [g l ipdE+ f fc/ipdE, 

J J u>0J J u<QJ 

where Wq = (po, mo, no, E 0 ) T . Since g l and g r have been obtained earlier, the above moments can be 
evaluated explicitly. Therefore, the conservative variables po,mo,no, and E 0 at the cell interface can be 
obtained, from which g 0 is uniquely determined. For example, Ao in go can be found from 

Ao = {K + 2)po/( / L(E 0 — 2^ m ° 2 n o)/Po ))• 

Then, d l and a r of g in Eq.(2.9) can be obtained through the relation of 


(2.15) 


and 


(2.16) 


/ar\ 

Wj+ijxj+i) - Wo _ - 0 d r 2 
PoAx+ a0 a r 3 

\d\) 


= M l 0 a 


0 ’ 


WQ-WjjXj) 

p 0 Ax- 


(d\ \ 
d 2 
d l 3 

\4 / 


= M a0<* 


l 

0 > 


where the matrix M® 0 = f go4’a4’0dE/ p 0 is known. Therefore, (o[, a 2 , a 3 , a 4 ) T and (a[,a l 2 ,a l 3 ,a l 4 ) T can be 
found following the procedure in Appendix B. 
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Up to this point, we have determined all parameters in the initial gas distribution function f 0 and the 
equilibrium state g at the beginning of each time step t = 0. After substituting Eq.(2.7) and Eq.(2.9) into 
Eq.(2.5), the gas distribution function / at a cell interface can be expressed as 

f(x j+1/2 ,t,u,v,0 = (l-e~ t/T )g 0 

+ (r(— 1 + e -t / r ) + fe -t / r j (a*H[ti] + d r (l - H[w])) ug 0 
+r(t/r - 1 + e~ t/T )Ag 0 

+e _i / T ((1 - u(t + T)a l )H[u]g l + (1 - u(t + r)o r )(l - H[u])<7 r ) 

(2.17) +e-‘/ T (-rA l E[u]g l - rA r ( 1 - U[u])g r ) . 

The only unknown left in the above expression is A. Since both / (Eq.(2.17)) and g (Eq.(2.9)) contain 
A, the integration of the conservation constraint Eq.(2.3) at Xj +1 / 2 over the whole time step At gives 

J J \g - f)i>dtdZ = 0, 

which goes to 

MlpAg = J [7i0o + 72 u (o'H[w] + o r (l - H[u])) g 0 

+ 73 (K[u]g l + (1 - R[u])g T ) 

+ 74 u {a l E[u\g l + a r (l - H[u])/) 

(2.18) + 75 ((o‘« + A l )E[u]g l + ( a r u + A r )(l - H[u])g r )] i/> a dE, 
where 


70 = At - r( 1 - e~^\ 

71 = -(l-e _At/r )/7o, 

72 = (-At + 2r(l - e- At / T ) - Me~ At > T ) / 7 o, 

73 = (1 - e~ At/T )/7o, 

74 = (A te~ At/r - r( 1 - e _A</r )) / 70 , 

75 = t{ 1 - e _At/T )/7 0 . 


Since all moments of the Maxwellian on the right hand side of Eq.(2.18) can be evaluated using Appendix 
A, Eq.(2.18) can be solved to find (.4, , ,4 2 , A3, At ) T using Appendix B again. 

Finally, the time-dependent numerical fluxes in the ^-direction across the cell interface can be computed 
as 


(2.19) 


(? P \ 

Fn 

\FeJ j+1/2 


/• 


= / U 


1 

u 

V 


\ 


f(x j+1/2 ,t,u,v,0dZ, 


\ |(u 2 + v 2 + £ 2 ) / 


where f(Xj + i/ 2 ,t,u,v,£) is given in Eq.(2.17). By integrating the above equation to the whole time step, 
we can get the total mass, momentum and energy transport. 
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2.4. Analysis. In this section, we are going to analyze the BGK scheme presented in the last section. 
Many issues related to the kinetic limits, collision time, Prandtl number fix, boundary condition, and kinetic 
model, will be addressed. 

2.4.1. Navier-Stokes Solver. In order to verify that Eq.(2.17) corresponds to a Navier-Stokes solu- 
tion, let’s consider the following limiting case. Eq.(2.17) gives explicitly the time-dependent gas distribution 
function / at the cell interface. In a well resolved flow region, such as in a resolved shock layer, the recon- 
structed conservative variables in Fig. (4.1) will become approximately a straight line. In such a case, the 
distribution function / 0 has g l = g r and a 1 = a r . Consequently, Eq.(2.14) gives g 0 = g l = g r , and Eq.(2.15) 
and (2.16) reduce to a 1 =a r = a 1 = a r . As a result, A determined in Eq.(2.18) is exactly equal to A 1 and A r 
in Eq.(2.13). Therefore, without any further assumption, the gas distribution function / at a cell interface 
becomes 

(2.20) / = go [l — r(ua + A) 4- t A\ , 

where —r(ud + A) go is exactly the nonequilibrium state in the Chapman-Enskog expansion of the BGK 
model [37], and go At is the time evolution part of the gas distribution function. The equation (2.20) is 
the equation we used for the low Mach number viscous flow calculations [29], where the accuracy of the 
above formulation is well established. Note that in deriving Eq.(2.20), we have not required the assumption 
r < At, which has been used previously [37]. The only requirement here is that the dissipative region is well 
resolved, such as the case with 5-10 grid points in the shock or boundary layers. In the under-resolved 
region, the BGK scheme will present a viscous solution for the discontinuous initial data. 

In the paper by Chae, Kim, and Rho [4], they basically interpolated goua as the nonequilibrium state of 
the Chapman-Enskog expansion. Actually, the correct form should be go{ua + A), and only this one could 
satisfy the requirement f go(au + A)tpdE = 0. 

2.4.2. Collision time. In a well resolved dissipative region, such as the cell size A.t is smaller than 
the dissipative length scale determined by the physical viscosity, the collision time r in the BGK scheme can 
be naturally determined by the physical relation 


t = p/p, 

where p is the dynamical viscosity coefficient and p is the pressure. The BGK model gives a fixed Prandtl 
number Pr = 1.0, which could only make one parameter correct, i.e., the viscosity or heat conduction. A 
numerical fix to make both coefficients correct will be addressed next in the Prandtl number fix part. For 
the viscosity coefficient, p can take any reasonable form in the determination of t . The simplest case is that 
p keeps a constant. In the shock-boundary interaction case, p will take the Sutherland’s law, 

P = P (— ) 3 / 2 r °° + S 

1 t + S ’ 

where T Xl and S are the temperatures with the values T rxi = 285 K and S = 110.4AT. 

Theoretically, the dissipative structure, such as the shock thickness, is solely determined by the physical 
viscosity. The structure should be independent of the cell size and time step used in a numerical scheme. 
However, even though the Navier-Stokes equations are accurately solved by the BGK method, if the cell size 
is not fine enough to resolve the wave structure, the physical one has to be replaced by a numerical one. For 
example, the physical shock thickness is replaced by the numerical cell size. In such a situation, we cannot 
solve the Navier-Stokes equations with the original physical viscosity. The effective viscosity in such a case 
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should be a combination of the physical and numerical ones. Different from many upwinding schemes, the 
BGK method cannot simply take the apology to admit that the implicit numerical viscosity is included in 
the under-resolved flow region. Since the BGK method is such an accurate Navier-Stokes flow solver, even 
in the under-resolved discontinuity region, the required additional numerical viscosity which is consistent 
with the numerical shock thickness has to be explicitly included. Since the jump in the flow variables at 
a cell interface, see Fig. (4.1), represents basically the under resolveness and appears automatically in high 
gradient flow region, the collision time r used in all simulations in this paper takes the following form, 

(2.21) r=^ + \ Pl/Xl ~ Pr/Xrl At, 

P | Pi / "-1 Pr/'^r\ 

where At is the CFL time step and the second part corresponds to the numerical viscosity. The second 
term on the right hand side in the above equation is related to the pressure jump at the cell interface in the 
reconstructed initial data. In the continuum flow region, this term will become very small. As shown in the 
test cases, this term neither poisons the boundary or shock layer calculations in the well resolved cases, nor 
reduces the shock capturing ability of the BGK scheme in the under-resolved region. 

The obvious advantage of the BGK method is that it solves a viscous governing equation with an explicit 
dissipative coefficient all the time, which avoids the ambiguity of implicit dissipation in many upwinding 
schemes due to the wave modeling in the Riemann solvers [31]. Even though the shock jump with a width of 
2 or 3 cell size can be captured nicely in the Godunov type schemes, the dissipation there for the construction 
of such a shock structure is solely coming from numerics. There is no reason to guarantee that the same 
numerical dissipative mechanism works in all physical situations [39]. The BGK scheme explicitly includes 
the physical and numerical ones into the algorithm. The fluid behavior in both smooth and discontinuous 
regions are described uniformly by the collisional BGK model with an adaptive collision time. The adaptation 
of collision time is necessary physically and numerically. 

Even the Navier-Stokes equations with an adaptive local viscous coefficient can be solved by the Godunov- 
type method, there are still difference between it and the BGK method. The BGK scheme gives a solution 
under the general initial condition (Fig.(4.1)) without separating the inviscid and viscous terms, it is difficult 
to design such an unsplitting time accurate Godunov method for the Navier-Stokes equations. Even for 
the same mass, momentum, and energy distributions, see Fig. (4.1) again, the kinetic scheme uses a non- 
equilibrium state /o to describe it. The macroscopic description could only see an equilibrium state initially. 

2.4.3. Limiting cases. The BGK scheme based on the BGK model is valid for the Navier-Stokes 
solution in both t < At and t > At region. If the collisionless Boltzmann equation f t + uf x = 0 is solved, 
the gas distribution function (2.17) goes to 

/ = fo(x - ut ) 

(2.22) = [l — r(ua l + A 1 ) — tua!] H [u\g l + [1 — r(ua r + A r ) — tua r ] (1 — H[u])(/ r . 

This is a 2nd-order time accurate scheme and is the limiting solution of Eq.(2.17) under the condition r 2> At. 
If the above distribution function is further simplified to the lst-order time accuracy, it becomes 

(2.23) / = [l — r{ua l + A 1 )] H[tt]r/ + [1 — r(ua r + .4 r )] (1 — H[u])< 7 r . 

The above distribution is ideally the same as the one used by Chou and Baganoff in their gas-kinetic Navier- 
Stokes solver [6]. In their approach, a direct implementation of the Chapman-Enskog distribution of the 
Boltzmann equation is used to split the flux. The non-equilibrium state in our case is solely consistent with 
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the BGK model. From the BGK scheme, we can dearly understand the limitation of the Chou-Baganoff’s 
Kinetic Flux Vector Splitting Navier-Stokes (KFVS NS) method. Because Eq.(2.23) is the limiting case of 
Eq.(2.17) with r 2> At, Eq.(2.23) is only valid for the Navier-Stokes solution under such a limiting condition. 
In other words, KFVS NS scheme approaches the Navier-Stokes solution accurately if the condition /i/p> At 
is satisfied. In p/p < At region, the KFVS NS scheme could behave badly for the Navier-Stokes calculation. 
In the later case, the free transport mechanism in the KFVS NS solver regards the time step At as the particle 
collision time, subsequently poisons the physical viscous solution. This artifact can be ignored only in the 
case At /i/p, where the physical viscous term is dominant. In order to get a more accurate understanding 
about the above analysis. In the following, we are going to qualitatively estimate the KFVS NS scheme in 
the shock and boundary layer simulations. Suppose we need N ~ 10 cells to resolve a NS shock structure or 
boundary layer. Since the shock thickness is proportional to the mean free path l s , in the shock layer case 
we need l s = N Ax. Then, the condition p/p At becomes 


— » At 


pcl s 

P 


» 


Ax 

\U\ + ~c 




M + 1’ 


where c is the sound speed, U is the macroscopic velocity, and M is the Mach number. For any shock 
wave with M > 1, the above relation can be satisfied. Therefore, the KFVS NS scheme could give an 
accurate NS shock structure [6]. We have also tested both Eq.(2.22) and (2.23) in the shock structure 
calculations in Case(l) of the numerical experiment section in the current paper. Both formulations give an 
accurate solution. However, if the boundary layer is resolved with the same number of grid points, we have 
\JvxjUoo = N Ax. In this case, the condition p/p > At goes to 


, r N 2 Ax 

M 

x 




1 

M + V 


where x is the distance between the point in the boundary measured and the leading edge, which is on the 
order of V 2 Ax. Therefore, for the subsonic boundary layer, such as M ~ 0.1, the above condition cannot be 
satisfied. The KFVS NS scheme cannot be an accurate NS solver in this case. Fig.(4.14) verifies the above 
analysis for both Eq.(2.22) and (2.23). Recently, it is interesting to observe that the KFVS NS method could 
give consistent results with DSMC simulation in the near continuum regime with the implementation of slip 
boundary condition [18]. The current BGK scheme can cover the similar cases. 

For the Euler solution, Eq.(2.23) can be further simplified. 


(2.24) 


/ = H [u]g l + (1 - H [u])g\ 


where the nonequilibrium state is totally removed. This is precisely the KFVS scheme for the compressible 
Euler equations [25, 26, 20, 23]. The above KFVS scheme has been well studied and applied to many 
physical and engineering problems. A earlier version of the above scheme is the beam scheme, where instead 
of Maxwellians the equilibrium states g l and g r are replaced by three Delta functions or particles [27]. As 
analyzed recently [30], the Steger- Warming method can be represented as a “beam scheme” too. But, due 
to their slight difference in the particle representation, such as the lack of internal energy in the second 
“particle” in the Steger- Warming method [30], it is less robust than the beam scheme. The relation between 
the beam scheme and the Lattice Boltzmann method is analyzed in [40]. 

With the above connection between the KFVS scheme and the Steger- Warming method, it is easy to 
understand the poor performance of many FVS schemes in the viscous boundary layer calculations [28, 33, 34]. 
Similar to the KFVS scheme, for the Navier-Stokes solution FVS methods also require p/p At. Due to 
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the symmetric lattice and diagonal transport, the Lattice Boltzmann Method (LBM) is very fortunate in this 
aspect. It does present an accurate NS solution in the incompressible limit [5, 12]. The reason for this is that 
with a symmetric lattice, the free particle transport mechanism from one node to another node could generate 
an artificial viscous term which is consistent with the Navier-Stokes term and its coefficient is proportional 
to — t. Therefore, in a fixed time step case the numerical dissipative term can be absorbed in the physical 
one [16]. As a result, the final viscosity coefficient in the Lattice BGK (LBGK) method is proportional to 
(r — At/2), where At = 1 is used there. There is no a precise analogue between the finite volume KFVS 
scheme and the Lattice Boltzmann method. Due to cell averaging, reconstruction process, and the non- 
isotropic transport, such as the lack of diagonal transport, the KFVS scheme has a much more complicated 
dissipative mechanism. But, the numerical viscosity coefficient v n can be still approximately estimated for 
a lst-order KFVS scheme using a simple shear flow model [19], which gives the same result as r = At in the 
LBGK method. The development of a multidimensional upwinding scheme will depend not only on the wave 
modeling, but more closely on the mesh construction. More precisely, it depends on whether a numerical 
mesh could preserve the isotropic and homogeneous properties of the fluid equations. CFD community 
usually has less experience in this aspect. There is something we can learn from the Lattice Boltzmann 
method, where the symmetry, invariants, etc., are the main concerns in their algorithm developments. In 
some sense, a triangular mesh has more symmetry and isotropic property than a rectangular one. 

2.4.4. Prandtl Number Fix. It is well known that the BGK scheme corresponds to unit Prandtl 
number. In order to change the above Prandtl number to any realistic value, many attempts have been 
proposed. The most well known one is the BGK-Ellipsoidal-Statistical (BGK-ES) collision operator [14], 
where the equilibrium state in the BGK model is replaced by an anisotropic Gaussian (without considering 
internal variables), 

G =7sm exp (- 1 i {u '- u,)r -' (u, - Ui> )' 

where pT = -p^RTI + (1 - A_)p<E> is a linear combination of the stress tensor /?<£ = / ( Ui - Ui)(uj — Uj)fdE 
and of the Maxwellian isotropic stress tensor pRTI. If we extend the current BGK scheme to the above 
BGK-ES model, considerable work has to be devoted to capture the time evolution of the above anisotropic 
stress tensor. 

As mentioned earlier, the BGK model itself can always make one coefficient correct, the viscosity or heat 
conduction. In the BGK method, we have obtained explicitly the time dependent gas distribution function 
/ at the cell interface Eq.(2.17). Therefore, the time-dependent heat flux can be evaluated precisely, 

(2.25) q=±J( u -U) ((« - U) 2 + {v- V) 2 + f ) fdE, 
where the average velocities U and V are defined by 

U = I ufdE/ 1 fdE , V= I vfdE/ j fdE. 

Then, the easiest way to fix the Prandtl number for the BGK scheme is to modify the energy flux by 
subtracting the above heat flux (2.25) and adding another amount with a correct Prandtl number, 

( 2 . 26 ) nr =fe + (^~ i )?, 

where TFe is the energy flux in Eq.(2.19). This fix can be equally applied to the BGK Discrete Velocity 
Model (DVM), where the discrete distribution function / is known [21]. 
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In a smooth flow region, the above Prandtl number fix can be further simplified. Since the gas distribution 
function in such a case reduces to / = go [l — r{au + A) + tA\ , the corresponding heat flux is 

q s = -t j g 0 (u - U 0 )(ip4 - U 0 ip 2 - V’oV’3 + * (I4 j + Vo)){du + A)dE 
= -T J go(u - U 0 )(ip4 ~ Uo4’2 - VoV’3 )(au + A)dZ 
(2.27) = —t J g 0 (au 2 ip 4 + Au ^’4 — Uo(au 3 + Au 2 ) — Vo(au 2 v + Auv))dE. 

So, in this case we can simply replace q in Eq.(2.26) by the above q s . This is the formula we are going to use 
in the Couette flow and the shock boundary interaction cases in the next section. There will not have much 
CPU time involved in the above Prandtl number fix, since all momentum in Eq.(2.27) have been obtained 
already in the evaluation of the original energy flux Fe . The above Prandtl number fix with the evaluation 
of q s is similar to the method proposed in [4]. The difference is that all terms related to A in Eq.(2.27) was 
ignored in [4], which will introduce errors in the unsteady flow calculation. As mentioned earlier, the correct 
nonequilibrium state in the BGK scheme should be proportional to ( ua + A) go instead of uago- 

The Prandtl number fix (2.26) is a post-processing correction, which is basically a numerical fix. But, 
to the Navier-Stokes order, the above fix is physically founded. Theoretically, the BGK-ES model is also 
a numerical fix, but it is on the dynamical level. Dynamically, keeping an anisotropic Gaussian for the 
equilibrium distribution function seem no any physical basis. The real physical weakness of the BGK model 
is that the collision time is independent of particle velocity, this fact is different from the phenomena with 
an anisotropic temperature distribution, where the temperature is directionly dependent. If only a correct 
Prandtl number is required, there is no reason to construct more and more complicated kinetic models. 

2.4.5. Boundary Condition. For the Navier-Stokes equations, the no-slip boundary condition is ob- 
tained by creating two ghost cells, where the velocities in the ghost cells are reversed from the velocities 
inside the computational domain, see Fig.(4.3). For the adiabatic wall condition, where there is no heat 
flow through the boundary, the mass and energy densities in the ghost cells should be symmetric around the 
boundary, 

P—i = Pi , E - 1 = Ei, 

P-i = Pi » E- 2 = E 2 , 

where —1 and —2 represent the 1st and 2nd ghost cells. Due to the above boundary condition, we can easily 
prove that the mass flux F p = 0 and the heat flux q = 0 at the wall for the BGK scheme. 

For the isothermal boundary condition, where the boundary keeps a fixed temperature, such as Ao = 
m/(2kT 0 ), in order to keep a 2nd-order accuracy of the scheme at the boundary we have to carefully derive 
the flow variables in the ghost cell. In the following, we only consider the case where the nonlinear limiter is 
not applied at the boundary cells. Therefore, we only need to construct the flow variables in the first ghost 
cell, see Fig. (4.3). In order to have the non-slip condition, we first have 

U-i = -lh , U-i = -V u 

which gives Uq = 0 and Vo = 0 at the wall. Since the temperature at x = 0 (location of the boundary) has 
a fixed value Ao, the slope of temperature in space at the boundary is 

( 0\ x _ Ai — Ao 

W°“ “pT’ 
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where Ai is the “temperature” in the first cell inside the computation domain. Therefore, the temperature 
in cell —1 becomes 

1 <9A 

A -' - Ao - 2 Ax( ^ )o - 

Since the energy density in cell —1 is determined by 
(2.28) = + 

the only unknown in the above formulation is p_ i. In order to determine p- 1 , we need to use the condition 
that there is no net mass flux transport across the boundary. Since the mass transport in a time step At 
can be expressed as (no limiter involved), 


which gives 


Am = 0 = J j u [l — r(au + A) + tA\ g 0 dtdE 

= J J u(l + t A) godtdE 


au 2 godr 


= At(p 0 U 0 ) - ^ (At) 2 J 


1A 0 = ±A 0 . 

po ox A 0 ox 


The discretized form of the above equation is 
(2.29) 


pi -p_i 1 d\ 


Pi + P—i 2Aq dx 


from which p_i can be derived. The above isothermal boundary boundary condition will be used in the 
Couette flow simulation. For a moving isothermal boundary, such as that with a velocity Vo, the only change 
from the above isothermal boundary condition is that V-velocity in the cell —1 is replaced by 

1 dV 

V ~' = 14 - 2 A *< & »»' 

where ( dV/dx) 0 = 2(\\ - V 0 )/Ax. 

Another important observation from the kinetic scheme is that it can introduce slip condition easily 
through the use of appropriate flux boundary condition. The kinetic scheme can be matched with the 
DSMC method in the near continuum regime only after the implementation of the slip condition at the solid 
boundaries. The basic formulation of kinetic slip boundary is based on the fact that with the introduction 
of gas distribution function, we can explicitly evaluate the amount of particles hitting the boundary, then 
according to the accommodation coefficients for the momentum and energy, and the temperature at the wall, 
we can re-emit the same amount of particles with a pre-described distribution function. As a result, the 
appearance of slip at the boundary is obtained naturally and is consistent with the DSMC type boundary 
condition in the near continuum regime due to their common kinetic considerations. More discussion and 
the specific application of slip boundary can be found in many kinetic books and papers [22, 17, 3, 2, 6, 18]. 
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2.5. Artificial viscosity — Godunov — BGK method. In the CFD algorithm development, the 
two classical pioneering papers for the shock capturing schemes are by von Neumann and Richtmyer [35] 
and by Godunov [8]. Since any physical solution has to be described in the discretized space and time, the 
limitation of cell size and time step has to be considered, von Neumann and Richtmyer realized that the 
numerical shock thickness needs to be compatible with the cell size. So, the central idea in [35] is that a 
viscous governing equation with an enhanced viscosity coefficient has to be solved numerically. 

The success in Godunov method is that it introduces a discontinuity in the flow representation. In the 
under-resolved flow simulation, due to the large cell size, a discontinuity will appear naturally in the initial 
data. The implementation of a discontinuity is much more important than the introduction of the Riemann 
solver. The cell interface discontinuity gives a more realistic representation about the physical situation. The 
numerical dissipation involved in the discontinuity can hardly be recovered by a delicate viscosity coefficient 
[39]. But, for a second order scheme with high-order initial interpolation, the dissipation introduced in 
the discontinuity is much reduced, which is not enough to construct a numerical shock wave. From our 
experience, we have the following conjecture. If a Generalized Riemann Problem (GRP) is correctly solved 
for the Euler equations with the inclusion of initial slopes in the gas evolution stage, see Fig. (4.1), a 2nd-order 
(in both space and time) accurate scheme cannot properly capture the numerical shock waves. Even with 
the discontinuity at a cell interface, additional numerical dissipation is still needed. As a special case, the 
Lax-Wendroff scheme is actually a generalized Riemann solver under the continuous initial condition. 

The methodology of the BGK scheme is in somehow to combine the two important issues risen in the 
above two methods, (i) a viscous governing equation with an enhanced viscosity coefficient (2.21) is solved, 
and (ii) follows the time evolution of the flow distribution from a discontinuious initial data. Both factors 
are important for the development of a robust scheme for the fluid simulation. In the smooth region, the 
additional numerical viscosity and the discontinuity at the cell interface disappear. The BGK scheme goes 
back to the traditional Lax-Wendroff type central schemes for the NS equations. 


3. Numerical Experiments. The current scheme has been applied to several test cases ranging from 
simple Couette flow to the complicated shock-boundary interaction case. Unless otherwise stated, in all 
numerical examples reported, the particle collision time is given by Eq.(2.21), 7 = 1.4. The time step At in 
all calculations are determined by CFL number equal to 0.7. All steady state solutions are obtained from 
the time accurate BGK solver with a long time integration. 

Case(l) Couette Flow with a Temperature Gradient 

Couette flow with a temperature gradient provides a good test for the BGK scheme to describe the 
viscous heat conducting flow. With the bottom wall fixed, the top boundary is moving at a speed U. The 
temperatures at the bottom and top are fixed with values T 0 and Ti . The analytic steady state temperature 
distribution is 


(3.1) 


T - To = V_ PrEc y _ y_ 

Ti-T 0 H 2 H' 


where H is the height of the channel, Pr is the Prandtl number, Ec is the Eckert number Ec = U 2 /C p (Ti —To), 
and C p is the specific heat ratio at constant pressure. 

We have set up the simulation as a ID problem in the x-direction. There are 20 grid points used in 
this direction from 0 to 1 with H = 1.0. The moving velocity at the right boundary in the y-direction is 
Vi = 1.0. The initial density and Mach number of the gas inside the channel are 1.0 and 0.1 respectively. 
The isothermal no-slip boundary conditions are implemented at both ends. We have tested the current BGK 
scheme with a wide range of parameters, (i). specific heat ratio 7 = 7/5, 5/3, and 2.0, (ii). different Prandtl 
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number Pr = 0.5,0.72,1.0,1.5, and 2.0, (iii). different collision time r, which ranges from 0.01 At to 1.5At. 
With the variations of parameters, all simulation results fit the exact solutions very well. Fig. (4.4) and (4.5) 
present the solutions in a few cases with different Prandtl numbers and Eckert numbers. From these figures, 
we see that the Prandtl number fix does modify the heat conduction term correctly. 

Case(2) Navier-Stokes Shock Structure 

This is also a ID case, which is mainly to test the performance of the BGK scheme from the shock 
capturing to shock structure calculation. The initial condition is that a stationary shock with Mach number 
M = 1.5 is located at x = 0. The viscosity coefficient for the Navier-Stokes equations takes a value v = 
0.00025, which corresponds approximately to a shock thickness, l s ~ 1/300. We have tested 6 cases with 
different cell size, which ranges from the under-resolved case with Ax = 1/100 to the well resolved case with 
Ax = 1/3200.0. The simulation results for density, velocity, and temperature distributions are shown in 
Fig.(4.6)-(4.8). From the coarse mesh to fine mesh cases, the shock structure gradually appears and it does 
converge to the exact Navier-Stokes solution. In the coarse mesh case, the BGK scheme could capture the 
shock jump crisply without any oscillation. Since the physical collision time r is determined by Eq.(2.21), 
with the change of cell size, the maximum ratio of t / A t ranges approximately from 0.3 to 10 through these 
cases. 

We have also tested the schemes based on Eq.(2.22) and (2.23) in the above shock structure calculation. 
Both methods give accurate Navier-Stokes solutions in this case, which is consistent with the observation in 
[ 6 ], 

Case(3) Mach 3 Step Problem 

The 2D Mach 3 step problem was first proposed by Woodward and Colella [36]. The computation is 
carried out on a uniform mesh with 120 x 40 cells, and the cell size used is Ax = Ay = 1/40. In order to test 
the viscous effect on the flow structure, we have used different Reynolds number Re = UL/v = 10 5 , 10 3 , 50, 
where the length scale is L = 1.0 and the upstream velocity is U = 3.0. The adiabatic slip Euler boundary 
condition is imposed at the boundaries in order to avoid the formation of viscous boundary layer. The density 
and pressure distributions at different Reynolds number are presented in Fig.(4.9)-(4.11). From these figures, 
we can clearly observe the effect of viscosity coefficient on the flow structure, such as the shear layer and 
the shock wave structure. Especially, in the case with Re = 50, the shock structure is well resolved. It is 
interesting to compare the BGK solution and the DSMC simulations in small Knudsen number regime, such 
as the case Fig. (14.29) in [2]. Due to the inclusion of non-equilibrium state and the easy implementation 
of kinetic slip boundary condition, the BGK scheme does provide a potential method to connect the Euler 
solution with the rarefied solution through the Navier-Stokes at near continuum regime. 

Case(4) Laminar Boundary Layer Case 

A laminar boundary layer with Mach number M = 0.15 and Re = 10 5 is tested over a flat plate. A 
rectangular mesh with 120 x 30 grid points is used and the mesh distribution in shown in Fig. (4.12). The 
U velocity contours at the steady state are shown in Fig. (4. 12). The U and V velocity distributions at the 
locations x = 6.438 and x = 34.469 are plotted in Fig.(4.13), where the solid lines are the exact Blasius 
solutions in the x and y directions. Due to the rectangular mesh, the number of grid points in the boundary 
layer at different locations are different. In both locations, the numerical solutions fit the exact solution very 
well. Due to the high Reynolds number in this case, the physical collision time r determined by the viscosity 
coefficient is much smaller than the time step At, i.e., r <C At. For this case, the previous BGK schemes 
could also capture the boundary solution correctly [37]. The results presented in Fig. (12) of [4] about the 
BGK scheme are due to the mis-use of the viscous term in the Chapman-Enskog expansion as stated earlier. 
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When Eq.(2.22) and (2.23) are used for the flux calculation in the above boundary layer case, the 
simulation results are shown in Fig. (4. 14). For the time accurate collisionless Boltzmann solution (2.22), the 
effective viscous coefficient is approximately proportional to ( Pphys + Afp/2), where p p h ys is the physical 
viscosity coefficient and At/2 is comming from the free transport mechanism. For the lst-order time accurate 
and 2nd-order space accurate scheme (2.23), the viscosity coefficient of the scheme is roughly proportional 
to (Uphys ~ Atp/2), where the forward Euler time steeping introduces an antidiffusive term (~ —A tp). 
Therefore, the solution in Fig. (4.14a) is more diffusive, and the solution in Fig. (4.14b) is less diffusive 
than the physical solution determined by p p h y s- Due to the stretched mesh, initial data reconstruction, 
and directional splitting, the motion of the numerical fluid from Eq.(2.22) and (2.23) has a complicated 
dissipative nature. The above effective viscosity estimates are only from the physical intuition. As we can 
see from Fig. (4. 14), the similarity solution is even lost. This means that the effective dissipative coefficient 
depends on local mesh size. 

Case(5) Shock Boundary Layer Interaction 

This test is about the interaction of an oblique shock at an angle 32.6° with a boundary layer. The Mach 
number of the shock wave is M = 2.0 and the Reynolds number for the upstream flow is Re = 2.96 x 10 5 . The 
dynamical viscosity p used here is the Sutherland’s law, and the Prandtl number is equal to 0.72. A mesh 
similar to Fig. (4. 12) with 110 x 60 grid points are constructed. The skin friction and pressure distributions 
at the surface of the plate is shown in Fig. (4. 15), where the data * is the experimental data from [9]. The 
pressure contours in the whole computational domain is shown in Fig.(4.16). Due to the high Reynolds 
number, the shock structure is not well resolved in this case. So, in terms of the shock, it is only a shock 
capturing scheme. But, in terms of the boundary layer, it is a Navier-Stokes solver because the boundary 
layer is well resolved. In this test, the condition r < At is satisfied. The previous BGK schemes work as 
well. The numerical solution from the previous BGK scheme is shown in [38]. However, for the KFVS and 
FVS schemes, the numerical dissipation could easily poison the boundary layer solution, since their validity 
condition is p/p^> At. 

4. Conclusion. This paper extends the BGK scheme to include the non-equilibrium state as the initial 
gas distribution function. As a result, a consistent BGK scheme for the Navier-Stokes equations is developed, 
which is valid in both r < At and r > At regions. The new scheme works not only in the viscous shear 
or boundary layer cases, where r < At is satisfied, but also in the construction of a Navier-Stokes shock 
structure, where r At. The KFVS NS method and KFVS scheme are the limiting cases of the current 
BGK scheme in the case p/p At, and both schemes are applicable to simulate the Navier-Stokes equations 
only under such a limiting condition. Also, the kinetic boundary condition, Prandtl number fix, as well as 
the relation among different schemes, are discussed in this paper. 

Following previous papers [24, 42, 41], the present paper shows a progressive development of the BGK- 
type schemes. The comprehensive numerical results presented in this paper and the physical and numerical 
analysis about the scheme indicate the level of maturity achieved by the gas-kinetic BGK method. 

Appendix A 

Moments of the Maxwellian Distribution Function. In the gas-kinetic scheme, we need to evaluate 
moments of a Maxwellian distribution function with bounded and unbounded integration limits. Here, we 
list some general formulas. 

Firstly, the Maxwellian distribution for a 2D flow is 

g = p^-^e-Miu-ur+iv-vf+e), 
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(u n+2 ) = U(u n+1 ) + n ^- 1 (u n ) 

L A 

When the integral is from 0 to +oo as (...) >0 or from — oo to 0 as (...) <0 , the error function and the 
complementary error function, appear in the formulation. Thus, the moments for u n in the half space are, 

<«°)> 0 = \edc(-V\U) 

1 e~ xlj2 

(u) >0 = U(u°) > o + -^- 
(u n+2 )> o = U(u n+1 ) >0 + =±l(«-)>o. 

and, 

(«°)<o = ^erfc{V\U) 

1 e~ xlj2 

(u} <0 =U(u°} <0 --^=- 
(u n+2 )<o = U(u n+1 ) <0 + ^l(«")<o 

Same formulation can be obtained for ( v m ) by changing U to V in the above moments of (it”). 


18 



Appendix B 

Solution of Matrix Equation b = M a . In the gas-kinetic scheme, the solution of the following equations 
is used many times, 



where b and M are known. The matrix M is from the integration of a Maxwellian, i.e., M a p = f ip a ippgdE/ p, 
and has the form 

/ 1 U V B\ \ 

U U 2 + 1 /2A( UV B 2 

M = 

V UV V 2 + 1/2A; B 3 ’ 

\Bi B -2 B 3 Bj 

where 

= ^(U 2 + V 2 + (K + 2)/2A), 

B2 = ^(*7 3 + V' 2 b r + {K + 4)U/2\), 

B z = \(V 3 + U 2 V + (K + 4)V/2\) 

and 

B 4 = * ((17 2 + E 2 ) 2 + (A + 4)(f7 2 + V 2 )/A + (K 2 + 6 K + 8) /4A 2 ) . 

The solution of Eq.(4.1) is the following. Define 

R 4 =2b 4 -(U 2 + V 2 + I ^)b 1 , 

R-3 = t>3 — V bi , 

i?2 = b ‘2 — Ub \ , 

the solution is 

04 = ZT2 (i?4 “ 2t7i?2 “ 2VRs) ' 

&3 = 2 A-R3 — V R 4 , 

CI2 — 2A/?2 — 

and 

«i = 61 — Ua 2 — V «3 — -a 4 (U 2 + V" H — — ). 
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x j+i/2 


Fig. 4.1. The 
can be constructed. 


reconstructed initial conservative variables around a cell interface, from which the nonequilibrium state fo 
The BGK scheme is based on the the solution of the collisional BGK model with the above initial condition. 



Fig. 4.2. The spatial distribution of the initial state fo and the equilibrium state g at t = 0. The evaluation 
fo is given in Eq.(2.1f)-(2.16). The final gas distribution function f in Eq.(2.17) at the cell interface Xj + 1/2 is a 
combination of fo and g. 


of g from 
nonlinear 
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Fig. 4.3. Boundary condition, (a), adiabatic boundary condition, where the mass, energy densities are distributed 
symmetrically around the boundary and the velocities are reversed, (b). isothermal boundary condition, where the velocity 
vector is reversed in the ghost cell, but the mass and energy densities are derived in Eq.(2.28) and (2.29). 



Fig. 4.4. Temperature ratio (T — Tq)/(T\ — To) in Couette flow. The solid line is the analytic solution given by Eq.(S.l), 
and the circles are the numerical results from the BGK scheme. While the Eckert number is fixed to 40, the Prandtl number 
takes the values 2.5, 1.0, and 0.72. The collision time r used in this scheme is about O.lAt. 
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Fig. 4.5. Continuation of Fig. (4-4 )■ While the Prandtl number is fixed to 0.5, the Eckert number takes the values 40.0, 20.0, 
and 4.0. 




Fig. 4.6. Density distribution of a stationary shock wave with M = 1.5. The kinetic viscosity coefficient of the flow is 
v = 0.00025, which corresponds to a shock thickness l a ~ 1/300. The numerical solution (+ sign) is obtained from the BGK 
scheme with different cell sizes. The solid lines are the exact Navier-Stokes solution. The cell sizes used are (a) 1/100, (b) 
1/200, (c) 1/400, (d) 1/800, (e) 1/1600, (f) 1/3200. 
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Fig. 4.7. Continuation of Fig. ( 4 . 6). Velocity distributions. 



Fig. 4.8. Continuation of Fig. (4-6). Temperature pj p distributions. 


















1 

0.8 
0.6 
0.4 
0.2 

1 

0.8 
0 6 
0.4 
0.2 

0.5 1 1.5 2 2.5 

Pressure 

Fig. 4.9. Density and pressure contours in Mach 3 step problem on a mesh with 120 X 40 grid points. The Reynolds 
number used in this case is Re= 10 5 w.r.t. the upstream velocity U = 3.0 and the channel height L = 1.0. The solution is very 
close to the solution of the Euler solvers. No special treatment is used around the step corner. 
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Fig. 4.12. (upper) Numerical mesh with 120 X 30 grid points for the boundary layer calculation, (lower) U velocity 
contours at iJe=10 5 . 
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